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The hydration free energies of ions exhibit an approximately quadratic dependence on the ionic 
charge, as predicted by the Born model. We analyze this behavior using second-order perturbation 
theory. This provides effective methods to calculating free energies from equilibrium computer 
simulations. The average and the fluctuation of the electrostatic potential at charge sites appear 
as the first coefficients in a Taylor expansion of the free energy of charging. Combining the data 
from different charge states {e.g., charged and uncharged) allows calculation of free-energy profiles 
as a function of the ionic charge. The first two Taylor coefficients of the free-energy profiles can 
be computed accurately from equilibrium simulations; but they are affected by a strong system-size 
dependence. We apply corrections for these finite-size effects by using Ewald lattice summation 
and adding the self-interactions consistently. An analogous procedure is used for reaction-field 
electrostatics. Results are presented for a model ion with methane-like Lennard- Jones parameters 
in SPC water. We find two very closely quadratic regimes with different parameters for positive and 
negative ions. We also studied the hydration free energy of potassium, calcium, fluoride, chloride, 
and bromide ions. We find negative ions to be solvated more strongly (as measured by hydration 
free energies) compared to positive ions of equal size, in agreement with experimental data. We 
ascribe this preference of negative ions to their strong interactions with water hydrogens, which can 
penetrate the ionic van der Waals shell without direct energetic penalty in the models used. In 
addition, we consistently find a positive electrostatic potential at the center of uncharged Lennard- 
Jones particles in water, which also favors negative ions. Regarding the effects of a finite system size, 
we show that even using only 16 water molecules it is possible to calculate accurately the hydration 
free energy of sodium if self-interactions are considered. 



" Address for Correspondence: Theoretical Biology and Biophysics Group T-10, MS K710, Los Alamos National 
Laboratory, Los Alamos, New Mexico 87545 

FAX: (505) 665-3493; Phone: (505) 665-1923; E-mail: hummer@tlO.lanl.gov 

^ Center for Nonlinear Studies, Los Alamos National Laboratory, MS B258, Los Alamos, New Mexico 87545 

Theoretical Chemistry and Molecular Physics Group T-12, MS B268, Los Alamos National Laboratory, Los Alamos, 
New Mexico 87545 



1 



I. INTRODUCTION 



A quadratic dependence on the ionic charge of the electrostatic free energy of solvation of a simple ion in aqueous 
solution is about the simplest reasonable possibility for that behavior. The Born model predicts that quadratic 
dependence.cl Several conwpiiter simulation calculations have shown that it is approximately correct for the simgksi 
monovalent ions in water fla Theoretical simplifications have been advanced to take advantage of such behavior.B'ETEl 

If that quadratic behavior were correct with sufficient accuracy, it would indeed permit important simplifications 
of the difficult task of molecular calculations of solvation free energies owing to electrostatic interactions in complex 
solutions. The theoretical simplifications identified on that basis can be viewed either as perturbation theory through 
second order in the electrostatic interactions, or as a Gaussian modeling of certain thermal fluctuations of those 
interactions. Adopting either view, these methods would have wide applicability and great simplicity. The question 
of the accuracy of the quadratic dependence on charge of the free energy owing to electrostatic interactions deserves 
to be raised for its own sake and given a precise answer as general as possible. 

This quadratic behavior is not a universal truth and previous simulation calculations have given helpful informa- 
tion on the circumstances where this quadratic dependence can be expected to fail.EI However, previous simulation 
calculations are sufficiently disparate that a high precision answer to the question of the accuracy of second-order per- 
turbation theory for the free energy owing to electrostatic interactions is not available. The disparate character of the 
available simulation results is largely caused by a lack of uniformity with respect to the treatment of finite-system-size 
effects on electrostatic interactions in aqueous solutions. It is not atypical for a finite-system-size correction and the 
electrostatic solvation free energy to be of similar size. 

In contrast to the role of compiuter experiments in answering this question, laboratory experiments have been useful 
mostly for framing the question.l3 Ej The difficulty of using laboratory experiments for the present purpose resides in 
our inability to extract generally an electrostatic contribution from contributions of the other interactions present. 

Because of these points, this work calculates the free energy owing to electrostatic interactions of simple, spherical 
ions in water by Monte Carlo methods and gives particular attention to the methodological issue of correction for 
finite system size. The molecular models used are simple but they have been widely tested. Because the goal of this 
work is to address the question of quadratic dependence on charge of the electrostatic solvation free energy, these 
models are sufficiently realistic for the present purposes. However, we will compare our computed free energies with 
experimental results and thus provide information on how these models might be simply improved for prediction of 
electrostatic free energies. 

Before proceeding with the technical developments it is worthwhile to give some discussion of the idea for the present 
treatment of system-size effects on solvation free energies of ions. There is no generally valid recipe that allows a 
determination of the effects of a finite system size on the calculated physical quantity in computer simulations. 
What must generally be done is to analyze the observed size dependence empirically. If, as is the case for Coulomb 
interactions of long range, different procedures are available, then we should expect consistent thermodynamic limiting 
(N — > oo) results for different methods of treating the finite-size system. It is well understood that certain quantities 
involving integrals over tiie whole sample, such as the dipole-moment fiuctuations, depend intrinsically on exterior 
conditions or constraints JHI Those conditions must then be properly understood theoretically. 

For the present problem involving the interactions and associated thermodynamics of an ion immersed in a dielectric 
liquid, a reasonable view is the following: Treatment of electrostatic interactions in a truly periodic format, e.g., by 
Ewald procedures, is consistent with the periodic boundary conditions that are nearly inevitable for other reasons. In 
periodic boundary conditions the interactions at the longest range that must be taken seriously occur at an appreciable 
fraction of the distance to the surface of the simulation cell. For typical simulated system sizes, ionic interactions at 
that longest range are large. Treatment of electrostatic interactions in a truly periodic format thoroughly tempers those 
large interactions. But a mathematical price for true periodicity of electrostatic interactions is a "self-interaction" 
associated with interactions with images and a uniform neutralizing charge background. For neutral systems this 
self-interaction can be sometimes ignored. For non-neutral systems such as those studied here there may be practical 
advantages of consistency obtained for explicit consideration of the self-interaction. We will account for these self- 
interactions explicitly in the calculations below. 

This argument permits treatments of the ionic interactions other than Ewald summation. In fact, the work below 
tests a generalized reaction-field (GRF) method and also finds that consistent results can be obtained if self-interactions 
are treated on a similar basis. 



II. THEORETICAL METHODS 
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A. Calculation of the free energy of charging 



The various methods to compute free energies using computer simulations have been rpjiewed extensively .Ii3~uil We 
start here from the potential distribution theorem for the excess chemical potential /i^^,E3 



^'"(91) -m'"(<Zo) = -fcBrin(exp{-/3[u(gi) -u(go)]})„ 



(1) 



where qo and qi are the two charge states and f3 = \/k^T] (. . .)q denotes a thermal configuration-space average in 
the charge-state g; and u{q) is the configuration-dependent interaction energy of the ion in charge-state q with the 
solution. Apart from finite-size corrections to be discussed later, u{q) is given by q(f>{r), where 0(r) is the electrostatic 
potential at the chargeposition r. i_. 

We next analyze eq |l| utilizing a cumulant expansiorJlj with respect to /3, 
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where Au = u{qi) — u{qo). This defines the cumulants C„ of order n = 0, 1, 2 as 
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We can interpret eq || as a Taylor expansion in Aq — qi — qo if we set Au = Aqcf) + (qf — qg)^/2, where ^ accounts for 
finite-size effects as a "self-interaction" to be discussed further below, 



A/i- = Aq ((</>) 
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where A/x*^^ — ^'^^{qi) ~ /i'^^((jo). The mean and the fluctuation of the electrostatic potential at the charge site q 
(corrected for finite-size effects) yield the derivatives of the free energy with respect to Aq. The information about the 
derivatives can therefore be extracted from equilibrium simulations. In principle, higher-order cumulants could be used 
to obtain information about the other Taylor coefficients. However, as was observed by Smith and van Gunsteren,cl 
higher-order cumulants are increasingly difficult to extract from computer simulations of limited duration. 

Therefore, we will evaluate Ci and C2 at few discrete charge states and combine this information about the 
derivatives, either by constructing an interpolating polynomial or by using a fit to a polynomial expression (or any 
other functional form) for the free energy as a function of Aq. The fit minimizes the mean square deviation of the 
observed data with respect to the coefiicients {0^} of the fitting function Afi'^^{q; {0^}), 



Afi^-{q,;{ak})-AfilUq,) 



A/i-((7,;{a4)-A/iXfe) 



(5) 



where Ui and pi are the estimated errors (standard deviations) of the observed first and second derivatives /i^^^ and 
/io^g at charge-state qi. 



B. Long-range Coulomb interactions and finite-size effects 

To minimize finite-size effects on energetic properties of Coulombic systems, we adopt the following strategy:0 We 
use lattice summation for calculating the electrostatic interactions to account for the periodic boundary conditions 
employed in the computer simulations; and we consistently include the self-interactions arising from lattice summation. 
We point out that aside from formal consistency the numerical results can motivate this approach by demonstrating 
in a finite-size analysis that the deviations from the thermodynamic limit (N —^ 00) are small. 

The Coulomb energy of a periodically replicated system of charges qi at positions (i = 1, . . . , N) can be expressed 

as 
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U= ^ qiqj ipEW (rij) + ^ ^ qUew : (6) 

l<i<j<N l<i<N 

where Yij = Vj — + n, with the lattice vector n chosen such that r^j is a vector in the upilj-ceil. The effective, 
position-dependent potential V3Ew(r) is obtained by lattice summation using Ewald's method,li3ll3'E2l 

, , erfc(?7|r + n|) .r-^ 47r f k'^ ■, \ 

- ? |r + n| + g V + • - VV ' 

where V is the volume of the box, erfc is the complementary error function, and fc = |k|. The two lattice sums extend 
over real- and Fourier-space lattice-vectors n and k, respectively. |— ■ |— ■ 

The self term ^ew = limr^o[</5Ew(r) — 1/?'] is the Wigncr potcntial:c2rEHl Using Green's theorem and A(l/r) = 
— 47r^(r), we find 



^EW = lim 



r-+0 



¥'Ew(r) - - 
r 



1 /■ 1 

= - — lim / dv - Av3Ew(r) • (8) 



r|>£ 

The integration region is infinite and includes all background charge and lattice image charges, 

1 



A(^Ew(r) = -47r 



E 



(9) 



Eqs ^ and ^ establish that ^ew is the electrostatic potential in a Wigncr lattice at a charge site owing to the 
lattice images-aud. the neutralizing background. For Ewald summation in a cubic lattice the self term is ^ew = 
-2.837297/Lj2rEl where L is the length of the cube. 

It will be interesting to remember that ^ew can also be expressed in terms of quantities associated with the primitive 
simulation cell 



i-ESN = -y- lim / dr i A(^ew(i") - y- / <y5Ew(r) n • V f i 



(10) 



The first term on the right is explicitly the interaction with the background density in the primitive simulation cell. 
The second term on the right is an integral over the surface of the primitive simulation cell. It describes electrostatic 
interactions of the central unit charge with a dipolar surface distribution </3Ew(r) n, where n is the surface normal 
pointing outwards. 

Eq 1^ can also be used for a non-neutral system since charges are implicitly compensated by a homogeneous back- 
ground in the Ewald formulation. This results in an expression for the energy difference Am between two configurations 
with different charge-states go and q\ of an ion at position r, 

Am = A9(^Ew(r) + ^^EW (g?-go) • (H) 

In the following, we will use this expression containing a self-interaction which is quadratic in the charge to-calculate 
the free energy of charging; i.e., we assume that the self-interaction accoimts for the finite-size correctionsEj 

In our calculations, we will also use a generalized reaction field (GRF).E3Ej The GRF Coulomb interaction depends 
only on the distance r of the charges and has a cutoff-distance rc, 

<PGRF(r) = ip(r/rc) e(rc-r) + C . (12) 
r 

O is the Heaviside unit-step function; p(a;) is a screening polynomial: 

p(x) = (l-x)'*(l-t-8x/5 + 2a;V5) . (13) 

By analogy with Ewald summation, we define the self term for the GRF as the potential at the charge site, ^grf — 
limr^o[¥'GRF(r) — 1/?"]. The total energy of neutral systems, if defined as in eq|^, is independent of the constant C. 
Howevac_in non- neutral systems C affects the total energy. We define C in analogy with the Ewald potential, which 



satisfie£3 

/ dv ^Ew(r) = , (14) 

JV 
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such that the average potential in the cell vanishes. If we require the normalization condition eq|lj also for the GRF 
interaction, we obtain C = —TTr'^/5V. The GRF self term is ^grf = — 12/5rc + C. For Vc = L/2, the normalization 
condition eq ^ accounts for only a small additional correction yielding ^grf = — 24/5L — 7r/20L. j— , 

It is interesting to make a connection with the correction method proposed by Sloth and S0rensen.tZl These authors 
use the minimum-image Coulomb interaction. To eliminate the system-size dependence in their calculation of chemical 
potentials of restricted-primitive-model ions, they introduce a background neutralizing the test-particle charge. This 
is done by adding a constant ^i/r to the bare Coulomb potentialJ23 

6/. = -^//rl. (15) 

This corresponds to enforcing eq ^ and adding a self term ^^/r = l™r^o[</'(f) ~ ^/^] for the minimum-image 
interaction, ^i/^ is also precisely the first term on the right side of eq It accounts for a large correction since 
Ci/r « -2.38/L. 



III. COMPUTER SIMULATIONS 

We calculated the free energy of charging ions in water using Metropolis Monte Carlo simulations. 00 The systems 
comprise a single ion and N water molecules. For water we used the simple point charge (SPC) modelEil The ion- 
water interactions were described by Coulomb and Lennard- Jones (LJ) interactions. The Coulomb terms involve the 
partial charges of oxygens and hydrogens on SPC water. The LJ interactions act only between water oxygen and 
the ion. We studied the ieas Na+, K+, Ca^+, F~, Cl~, and Br~. The LJ parameters for these ions were those of 
Straatsma and BereadsenES We also studied the chaxging of a model ion Me with methane LJ parameters as given 
by Jorgensen et alS3 Lorentz-Berthelot mixing rulcsEj were applied to obtain LJ parameters between water and Me. 
The L J parameters are compiled in Table ||. 

The charge interactions in the simulations were calculated using Ewald lattice summation (eqs ^ and ^ and the 
generalized reaction-field potential (eqs ^ and |l^. In both cases, the real-space interactions were truncated on an 
atom basis using L/2 as cutoff and applying the periodic boundary conditions on an atom basis. For the Ewald 
Fourier-space calculation, a cutoff < 38(27r/L)^ was used resulting in 2 x 510 k vectors. To correct the background 
dielectric constant from infinity to erf = 65, a term 27rM^/(2eRF + 1)^^ was added to the potential energy (in both 
Ewald and GRF calculations), where M is the net dipole moment of the water molecules. The real-space damping 
factor was set to 77 = 5.6/L. Electrostatic potentials at the ion sites were calculated using tpEW and fORF- The 
potentials were calculated after each pass (one attempted move per particle) and stored for analysis. For each system 
100 000 passes were used for averaging. Random configurations or configurations of previous runs were used as 
initial structures and always extensively equilibrated. The temperature was 298 K. The total number density was 
p = 33.33 nm~'^ in all simulations. Cubic boxes were used as simulation cells with edges L — [{N + l)/p]^/^. The 
Monte Carlo move widths were chosen so that an approximate acceptance ratio of 0.5 was obtained. 

In addition, thermodynamic integration (TI) was used to calculate directly the free energy of charging. Within 
100 000 Monte Carlo passes, the charge of the ion was linearly changed from to its full magnitude (±e, 2e, where e 
is the elementary charge). The free-energy changes were then calculated as 

n 

p^%q{) - p^^iqo) = {qi - qoh'^ ^^l' ' '?o)/2 , (16) 

where the sum extends over n = 100 000 Monte Carlo passes and the last term is a finite-size correction. Eq |l^ 
approximates the exact expression 

TI was also performed using the reverse path, i.e., decreasing the charge to 0. 

We also performed Monte Carlo simulations of ion-water clusters comprising one ion and SPC water molecules 
(4 < iV < 256). The starting structure was a random configuration with bulk density of water in a cubic box around 
the ion. The cluster was equilibrated for at least 50 000 passes (with an acceptance rate of about 0.5) and then 
averaged over 50 000 passes at 298 K. We used the bare Coulomb interaction l/r and did not apply a distance cutoff. 
No periodic boundary conditions were employed in the cluster simulations. 
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IV. RESULTS AND DISCUSSION 



A. Charging of a methane-like Lennard-Jones particle 



The free energy of charging a methane-hke LJ particle in SPC water was determined from a series of simulations 
with N = 128 and 256 water molecules and with Ewald and GRF charge treatment. A range of charges from — e to 
e was covered in steps of 0.25e (A'' = 128) and 0.5e (TV = 256). The results for the mean m and the fluctuation / of 
the potential at the ion site (with and without finite-size correction) are compiled in Table ||. In the calculations, the 
potential (p at the ion site (r=0) is defined as 



N 



= I]I]%„'P(r»J , (18) 

where the double sum extends over all water oxygen and hydrogen sites; ip is either (^ew or '/'grf- The mean m and 
the fluctuation / are calculated from 100000 Monte Carlo passes as 

m — e{<f)) , (19a) 
/ = /3e2((0_(0))2) . (19b) 

The corrected values for mean and fluctuation are defined as nic = m + qeS^ and fc — f ~ e^^- The Taylor expansion 
of the free energy of charging around a charge-state q assumes the following form: 

From the results of Table || we see that the finite-size corrections are of similar magnitude as the uncorrected results m 
and /. The uncorrected results of the different methods and system sizes are widely spread. If however the finite-size 
corrections are applied, we obtain consistent results for all methods and system sizes over the range of ion charges 
considered. With estimated errors (one standard deviation, as calculated from block averages) of 4.0 and 30 kJ mol~^ 
for nic and fc, we find data of different methods within two standard deviations from each other. The only exception 
is the fluctuation fc for q = — e, where the two extreme values (Ewald, N — 128 and 256) differ by about three 
standard deviations. In the following, we will restrict the discussion to the corrected values ruc and fc- 

Figure |l] shows the probability distribution P{ecj)) of the electrostatic energy ecj) for an ionic charge of q — e, 
as calculated from histograms. The P{e(t)) curves follow closely Gaussian distributions with mean and variance 
calculated from the </> data. This reflects the approximate validity of second-order perturbation theory in the ionic 
charge. However, the P{e(t>) curves for Ewald summation with N = 128 and 256, as well as GRF with N — 256 water 
molecules differ widely, both in the peak position and in the width. To illustrate the importance of the flnite-size 
correction, we included in Figure ^ the Gaussian distributions corresponding to the corrected values rUc and fc for 
mean and variance. The application of the finite-size corrections brings the three curves to very close agreement, 
yielding results that are approximately independent of system size and treatment of electrostatic interactions. 

To further illustrate the importance of the finite-size correction, we calculated (0) from the pair correlations of the 
Ewald-summation simulation with N — 256 water molecules as 

{<I)){R) = 47rpH20 / dr (^(r) \qogio(r) + 2gH.giH(r)] . (21) 

PH2O is the water density, qo and are the oxygen and hydrogen charge, and g\Q and gin are the ion-oxygen and 
ion-hydrogen pair correlation functions. Figure g shows the results for the charge-state g = — e of the ion Me as a 
function of the integration cutoff R for the bare Coulomb potential </?(r) = 1/r and (^grf with Vc — L/2. In both cases 
we included the flnite-size correction as a constant. The integration of the 1/r interaction extended into the corners 
of the cube, using the correct weights. As a reference, the Ewald result is shown as a straight line. All three methods 
converge to within about 1 kJ mol^^, which has to be compared with the estimated statistical error of 4 kJ mol^^ of 
the data. The integrated 1/r interaction shows strong oscillations and only in the corners of the cube does it approach 
its flnal value. The GRF interaction on the other hand contains a large self term and within two oscillations reaches 
its limiting value. 

This illustrates an important point regarding the correction of finite-size effects in the calculation of charge-related 
quantities. We achieve agreement between different methods of treating Coulomb interactions (Ewald summation, 
reaction field, bare Coulomb interaction) if we (i) normalize according to eq [ij and (ii) add a self term ^ = 
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limr^o[ip{r) — 1/r] to the energy. Further demonstrations of the vaHdity of these finite-size corrections will be given 
in the discussion of the results for sodium and fluoride ions in SPC water. 

Figure || shows rric as a function of the charge. We observe two linear regimes with different characteristics for g < 
and q > 0. Linear behavior of rric on the whole range of q would reflect validity of the second-order perturbation 
theory. It would imply Gaussian statistics of (jj and, correspondingly, that the coefficients in the Taylor expansion 
of order three and higher vanish. However, since we observe a transition in the linear behavior between charges of 
— 0.25e and 0, the statistics are only approximately Gaussian. We note that from the data of 100 000 passes it 
proved impossible to extract reliable information about the Taylor coefficients (cumulants) of order three and higher. 
The second Taylor coefficient fc can however be extracted accurately. Figure^ shows fc as a function of q/e. Included 
in Figure ^ as lines are the values of fc estimated from the linear fits of rric for q < and q > 0. 

We have fitted the rUc and fc data by a model with two Gaussian regimes. Included in Figures ^ and ^ is a fit 
of the whole set of derivative data (38 data points) to 

H'^'^iq) - /i"^(0) {a+q + b+q^)[l + tanh(c + dq)]/2 + [a^q + h^q^)[\ - tanh(c + dq)]/2 , (22) 

where is defined as in eq ^ with parameters a+, &+, a_, &„, c, and d. This model can nicely reproduce the data. We 
find a transition at q = c/d ~ — 0.2e between the two regimes of approximately Gaussian behavior with a quadratic q 
dependence. We ascribe this transition to differences in the structural organization of water molecules near negatively 
and positively charged ions. A possible explanation for the observed behavior is that for positive ions, the oxygen 
atom of water is pointing towards the LJ particle. The strongly repulsive forces of the r~^^ interaction prevent large 
fluctuations of because of the restricted oxygen motions. The hydrogens are pointing away so that rearranging 
them has only a comparably small effect on 0. For negative ions, the structures with one of the hydrogens pointing 
towards the ion will dominate. Because of the symmetry between the water hydrogens and the finite life time of the 
hydration shell, transitions will occur which could explain the larger fluctuations in the negative charge range. 

Similarly, a transition to a different Gaussian behavior for highly-charged positive ions was observed by Jayaram et 
al.a These authors studied the free energy of charging of a sodium ion in the charge range to 3e. When increasing 
the ion charge, a transition occurs to a more weakly decreasing-quadratic free-energy regime at a charge of about 
l.le. This transition has also been discussed by Figueirido et al£3 r-. 

We also find a nonvanishing potential at the methane site even at zero charge.B In a dipolar solvent, ((/))q=o is zero 
because of charge-reversal symmetry. However, if higher multipole moments are present on the solvent molecule, this 
symmetry is not conserved. The asymmetry of the charge distribution on the water molecule gives rise to a positive 
potential for q — 0; this is primarily caused by the hydrogens penetrating the LJ sphere of the methane particle, since 
they do not have a protecting repulsive shell in the model used. As a consequence, there is a small charge region in 
which increasing the charge Mosts free energy. A positive potential at the center of an uncharged particle was also 
observed by Rick and Berne .E3 

As a consequence of both the positive potential at zero charge and the larger potential fluctuations for negative 



ions, negative ions are more stably solvated compared to positive ions. Table [H compiles the free energies of charging 
as calculated from fitted polynomials p„ of degree n to the derivative data TOc and fc- Except for the simple Gaussian 
model p2, different fitting functions give consistent results for the free energies of charging. For ions with charge e 
and — e we find A/i'^^ = —250 and —431 kJ mol~^. Interpreted within a Born model for the free energy,tl i.e., 

A/ig^orn = -(1 - l/e)9V2i? , (23) 

we obtain Born radii i?+ = 0.27 nm and i?_ = 0.16 nm. (A value of e = 80 is used for the dielectric constant, but this 
hardly affects the results). The difference between i?+ and i?_ is somewhat smaller if we use the actual coefficients 
of the q^ term in the free energy expansion, as obtained from eq |2^, giving 0.23 and 0.16 nm for the Born radii of 
positive and negative ions. We emphasize the model character of the interaction potentials used in this study. A 
repulsive shell of the hydrogen atom might reduce the free energy difference between positive and negative ions. The 
favoring of negative ions however should persist. 

The lower free energy of negative ions compared to positive rions of equal size agrees with the experimental ob- 
servations. The hydration free-energy data compiled by Marcua23 for alkali metal and halide ions show a power-law 
dependence with respect to the ion radius. Using these fitted curves, we find differences of 150 and 240 kJ mol~^ for 
the solvation free energy between negative and positive ions of the size of potassium and sodium, respectively. The 
LJ particle Me studied here has a van der Waals radius between those of K+ and Na+. The calculated free energy 
required to go from — e to +e is 180 kJ mol^i. which is indeed bracketed by the experimental data. 

The revised Born model by Latimer et aZ.LJ also yields lower free energies for negative ions. For alkali and halide 
ions, it uses effective Born radii i? = rp + A, where rp is the Pauling radius and A is 0.085 and 0.010 nm for cations 
and anions. This smaller effective-radius correction for anions in eq results in considerably lower free energies of 
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negative ions compared to positive ions of equal size-in agreement with our calculations. The difference of the effective 
Born-radius correction as defined by Latimer et aLt3 is 0.075 nm, which agrees with what we find for the Me ion. 

The energetic differences in the hydration of positive and negative ions go along with differences in the structural 
organization of water molecules in the hydration shell. Figure |^ shows the ion-water pair correlation functions for 
different ionic charges. Going from g = to positive charges does not change the qualitative properties of the ion- 
oxygen and ion- hydrogen correlation functions gio and gia. An increase of the ionic charge results in a higher first 
peak. However, going from charge g = to negative charges affects strongly the structural organization of the first 
hydration shell. Already at g = — 0.5e, giYv shows the buildup of a second peak at about r = 0.2 nm distance. This 
peak reaches a value of almost 5 at g = — e, compared to gm essentially being zero in this distance region for charge 
g = 0. This strong interaction of the negatively charged ion with the hydrogens of water in turn affects the ion-oxygen 
correlation functions. Despite the negative charge of both the ion and oxygen site, gio has a first peak with a height 
of about 5 for g = — e compared to only 3 for q — e. The strong charge repulsion between water oxygens and the ion 
with g = — e is overcome by a large attraction caused by a water hydrogen pointing towards the ion and penetrating 
the ionic van der Waals shell without energetic penalty. 



B. Free energy of charging of the ions Na+, K+, Ca2+, F", CI", and Br" 

Using the L J parameters of Straatsma and Berendsenll (see Table ^ , we computed solvation free energies of ions 
representing Na+, K+, Ca^+j F~, CI", and Br". A gain, we emphasize the model character of this study. Its purpose 
is not to provide accurate theoretical values for the free energies but rather to characterize the theory. We can 
expect to obtain accurate values only after considerable improvement of the currently rather crude descriptions of the 
interaction potentials used here, aad. similarly in most other studies. Some of that work has indeed been guided by 
using free energies of hydration.ol2j However, controversies about certain technical aspects, ipriHiarily regarding the 
correct treatment of long-range interactions, need to be resolved to obtain conclusive results.E3cJ 

We extensively studied the solvation free energy of the sodium cation using the model described in section IH. 
Monte Carlo simulations using N — 128 water molecules were carried out for charges 0, 0.5e, and l.Oe to calculate the 
mean nic and the fluctuation fc of the electrostatic potential cj) at the ion site. As in the previous calculations, 100 000 



passes were used for averaging. The results are listed in Table IV. As for the uncharged methane, the potential at 
the uncharged sodium site is slightly positive. The decrease of mj, with increasing charge is stronger than linear and, 
correspondingly, the fluctuation fc increases slightly with the charge. This indicates that a simple Gaussian model 
using an expansion around the uncharged particle is of limited utility. 

We use the information about the derivatives to calculate the free energy of charging using polynomial fits. The 
results for the sodium ion using polynomials of degree 2, 4, and 6 are compiled in Table^ Also included in Table ^ 



are results obtained from TI, as described in section IH. TI was performed using Ewald summation and N = S, 



16, 32, 64, 128, and 256 water molecules as well as using the GRF Coulomb interaction and N = 32, 64, and 128 
water molecules. We observe excellent agreement of the free-energy data from polynomial fits and TI, except for the 
P2 fit which cannot fully account for the increasing potential fluctuations with increasing charge. The TI data of 
charging from to e and uncharging from e to show variations of about 5 kJ mol"^. Regarding the system-size 
dependence, Ewald summation gives accurate results even for as few as = 16 water molecules. The GRF shows a 
more pronounced system-size dependence with the = 64 data (cutoff Tc = 0.62 nm) being slightly too low. These 
results indicate that the free energy of charging is unexpectedly insensitive to the system size if the electrostatic 
interactions are treated appropriately. In particular, it is important to apply the correct finite-size corrections. For 
Ewald summation with N — 16, for instance, the finite-size correction accounts for about 60 percent of the free energy. 
Without the self terms the Ewald results for A^ = 256 and A^ = 16 differ by about 63 kJ mol"^; with the self terms 
included the difference is only 5 kJ mol~^. 



Table VI lists the results of polynomial fits of the free energy to the derivative data for the other ions studied (K"*" , 
Ca^"*", F~, Cl~, and Br"). Also included are results of TI calculations using Ewald summation and A^ = 128 water 
molecules. Except for the polynomial fit of degree 2, we obtain consistent results from the derivative data and TI. 
The p2 results are always somewhat too negative but this is more apparent for the negative ions. The two TI data 
per ion typically bracket the p4 and pe results for the free energy. 

Interestingly, there is no simple trend for the free energy of charging of monovalent cations with the ion size (as 
measured by a of the LJ interaction). The positive ions Na"*" and K+ as well as the negative ions F", Cl~, and Br~ 
show the expected decrease of A^*^^ with increasing a. However, only the negatively charged methane-like LJ particle 
Me~ fits into this ordering. The positively charged Me+ has a less negative A^*^^ than K"*", even though the van der 
Waals diameter a of K+ is considerably larger. However, the LJ interaction of the K+ ion is more shallow than that 
of Me"*" with the LJ e values differing by a factor of about 150. 



8 



We also calculated the excess chemical potential of inserting uncharged LJ particles iii_|SP|C, water of density 
p = 33.33 nm-3 at temperature T = 298 K. This was done using test-particle insertion.Ea&El 5000 SPC water 
configurations were used of a simulation run extending over 500 000 Monte Carlo passes. The simulation was performed 
using TV = 256 water molecules and GRF Coulomb interaction with a cutoff of Tc = 0.9 nm. We calculated (exp(— 
using 100 test particles per configuration, where u is the interaction energy of a LJ test particle with the water 
molecules. For the LJ interaction, a spherical cutoff distance of L/2 = 0.9865 nm was used. A cutoff correction for 
the term was applied, assuming homogeneous water density beyond the cutoff. The excess chemical potential is 
calculated as 

= -fcBTln(exp(-/3M)) . (24) 



Results are listed in Table VII. We find positive values for /i'^^ between 9 and 25 kJ mol^^, favoring the gaseous state. 
Adding /i*^^ to the free energy of charging, we obtain single-ion free energies of hydration. 

Experimental data_£or single-ion free energies OiLJiydration have been compiled by, for instance, Friedman and 
Krishnan.Eil Conway,E3 and most recently Marcus.Ej The first two references report values for the standard molar 
Gibbs free energy AG°, i.e., for a hypothetical transfer from a 1 atm gas state to a 1 mol/1 solution. Marcus lists 
values for AG* which is the Gibbs free energy of bringing an ion from an empty box into solution. The theoretical 
calculations determine the excess free energy of hydration, i.e., the transfer from an ideal gas of given density to 
solution of equivalent solute density. This process corresponds to that used by Marcus, so that AG* is the experimental 
equivalent of the theoretical free energy that we have referred to as /i^^ disregarding volume contributions. Because 
Marcus used AG* for the experimental free energies of hydration we will retain that notation here for those quantities. 
To convert from AG° to AG*, requires adjustment for the differences in standard states: we add to AG° the free 
energy of an ideal gas going from a pressure po corresponding to a density of 1 mol/1 to a pressure pi = 1 atm, 
which is A:BTln(pi/po), i.e., AG* = AG° — 7.92 kJ mol^^.Ea-Another correction accounts for differing values for the 
reference ion H+. We take the most recent value by Marcua^j AG*[H''"] = —1050 ± 6 kJ mol^^ and adjust the other 
values [-1098 (ref |l|) and -1074 ± 17 kJ mol"^ (ref ||)] accordingly. 

Results for the calculated free e nergy of ionic hydration /i*^^ = ii'^^{q = 0) + A/i'^^(0 q) and the experimental 



values AG* are compiled in Tabic VIIL For the calculated values we use those obtained from a fit of a sixth-order 



polynomial pg to the derivative data, as listed in Table VL The experimental data were adjusted as described above. 



The experimental data for cations show little variation between the three sources. However, the anion data vary by 
as much as 70 kJ_ixLol^^, with the Conway dataEd bracketed by the those of refs ^ and |5l|, but generally closer to the 
data of Marcus.L^I 

The calculated free energy data for cations do not show a clear trend. The results for Na+ and K+ are too low and 
too high by about 10 percent, respectively. The hydration free energy of Ca^+ is too high by about 15 percent. The 
anions on the other hand show a clear tendency with the magnitudes of the calculated free energies generally being 
too large. The relative errors are 26, 10, and 15 percent for F~, Cl~, and Br~, respectively, compared to the data 
of Marcus. These significantly too negative values of the hydration free energy of anions might be a consequence of 
the unprotected hydrogen atoms in the water-ion interaction model used. The positively charged hydrogen atom can 
penetrate the LJ shell of the ions without a direct energetic penalty. The interaction with the negative point charge 
at the center of the ion strongly binds the water molecule, resulting in a large enthalpic contribution to the free energy 
of hydration. But also effects of non-additive interactions might play a considerable rolea 



Also included in Table VIIl arc computer simulation results by Straatsma and BerendsenH These authors used 
thermodynamic integration in conjunction with isothermal-isobaric molecular dynamics simulations to compute hy- 
dration free energies of ions. The interaction potentials used here are identical with those of Straatsma and Berendsen, 
except for the treatment of the electrostatic interactions. We used Ewald summation; whereas Straatsrria-and Berend- 
sen used a spherical cutoff and a Born-type correction for finite-size effects. These authors (and otherfj) argue that 
the application of a Born-type correction is rather crude, approximating the solvent molecules beyond the cutoff by|-a 
dielectric continuum. Nevertheless, in the absence of a better alternative it has been widely adopted. Migliore et aZ.Ej 
calculated the free energy of ionic hydration based on a perturbation formula from Monte Carlo simu lations using 



MCY water and ah initio ion-water potentials. These authors also used a spherical cutoff. Table VIII includes the 
results of Migliore et al., who did not apply a finite-size correction. _ _ 

Qualitatively, our free-energy data agree with those of Straatsma and Berendseno and Migliore et a/j£j We observe 
the same ordering of the free energies with respect to ion size. The quantitative agreement is however poor. Our values 
for the cations Na+ and K+ are closer to the experimental data of Marcus. The cation free energies of Straatsma 
and Berendsen (with Born correction) are consistently more negative than those of our calculations. On the other 
hand, our anion free energies are significantly more negative than those of Straatsma and Berendsen as well as of 
Migliore et al. The results of Straatsma and Berendsen for Cl^ and Br^ are somewhat closer to the experimental data 
of Marcus when the Born correction is included. Without the correction they are significantly too high. The most 
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pronounced discrepancies between the anion data of Straatsma and BerendsenE3 and ours are those of the fluoride 
ion, with our fi'^^ values being lower by 83 kJ mol~^. This is somewhat surprising since Straatsma and Berendsen 
used the same parameters for the water-water and water-ion interactions. The difference can be a consequence of 
using different ensembles (NVT versus quasi- NPT); or, more likely, it is caused by the different treatment of the 
electrostatic interactions (Ewald versus spherical cutoff). 

The fluoride ion also shows the largest relative deviations from the experimental results. To further investigate 
these discrepancies, we have studied the energetics of clusters of different size formed by a single fluoride ion and 
water. We have performed Monte Carlo simulations using one ion that nucleates TV = 4, 8, 12, 16, 32, 64, 128, and 
256 water molecules at 298 K, as described in section [11. We calculated the interaction energy Ug of the fluoride ion 
with the SPC water molecules over 50 000 passes. Figure^ shows the differences Aw^ = Us — Us,ew with respect to 
the bulk simulation using N = 128 water molecules and Ewald summation as a function of N^^/^ . Aug can be fltted 
to a third-order polynomial in N^^^^ over the whole range of system sizes. Extrapolation to — > oo yields a limit 
for Aus close to zero. (However, the nontrivial dependence on N~^/^ limits the accuracy of the extrapolation.) The 
result obtained from Ewald summation, Ms,ew = — 1077±4 kJ mol"^, also agrees with the value Ug = —1075 kJ mol~^ 
obtained from integrating the pair correlation functions of the bulk simulation using (^(r) = 1/r in eq ^l], adding the 
LJ contributions, and applying the finite-size correction — e^^i/^- The integration shows that the LJ contributions are 
strongly repulsive (« 90 kJ mol"^) but compensated by large electrostatic interactions. 

The value for the solvation energy reported by Straatsma and Berendsen, Ug = —823 kJ mol^^, is however con- 
siderably smaller. The observed differences in Ug of about 150 kJ mol^^ agree in magnitude and sign with those of 
the free energies (83 kJ mol~^). If we truncate the integration of 1/r weighted with the pair correlation functions 
obtained from Ewald summation at i? = 0.9 nm (which is the cutoff Straatsma and Berendsen used) and do not apply 
a finite-size correction, we obtain a value of —867 kJ mol^^ in much closer agreement with Straatsma and Berendsen's. 
This indeed indicates that the treatment of the electrostatic interactions (Ewald summation versus spherical cutoff) 
is the major source of the discrepancy. 

Also included in Figure |^ are the results for the mean and the variance of the electrostatic potential at the ion 
site. Figure ^ shows differences with respect to the bulk value. The differences of the mean values A((^) closely follow 
the solvation-energy differences Awg and can also be fitted to a third-order polynomial in iV~^/'^. The differences 
of the fluctuation A(A0^) depend linearly on N"^/^ for N between 8 and 256. Both fitted curves extrapolate to 
approximately 0, indicating that the calculated bulk values are the correct limits for N oo. 

From the cluster-size dependence of the solvation energy and the mean and variance of the electrostatic potential, 
as well as the results for Me and Na+, we conclude that the use of periodic boundary conditions in conjunction 
with Ewald-summation (or reaction-field) electrostatics closely approximates the correct bulk behavior of the system; 
however, to get correct energetics it is important to include the self-interactions in the Coulomb energy. 



V. CONCLUSIONS 



We have shown that free energies can be accurately calculated from equilibrium simulations by extracting derivative 
information with respect to a coupling parameter. We have studied the free energy of charging ions in water, which 
accounts for most of the free energy of ionic solvation for typical ion sizes. The choice of the ionic charge as coupling 
parameter results in free-energy expressions involving cumulants of the electrostatic potential 4> at the charge sites. 
We find that the statistics of (j> are approximately Gaussian. This means that only the first and second moment of 
the distribution can be calculated accurately, with higher moments dominated by the poorly sampled tails. Corre- 
spondingly, only information about the first and second derivative of the free energy can be calculated accurately for 
any given charge state. The information for different charge states {e.g., uncharged and fully charged) can then be 
combined using interpolation or polynomial fitting. 

We have studied a methane-like Lennard- Jones particle in SPC water. We observe two almost Gaussian regimes 
separated by g = with different characteristics. Negatiiffi ions are more stably solvated compared to positive ions 
of equal size, in agreement with the experimental data.Ej The system shows further asymmetry, since the average 
electrostatic potential at the position of the uncharged particle is positive. This means that increasing the ion charge 
first costs energy. We relate these asymmetries of the energetics (lower free energy of negative ions, positive potential) 
to the structural asymmetry of the water molecule. The hydrogen atoms can penetrate the ionic van der Waals shell, 
whereas the oxygen atom is better protected. For the uncharged particle, this results in a net positive potential; and 
the point charge at the center of negative ions exerts strong electrostatic interactions with the tightly bound hydrogen 
of water. 

However, particularly for small anions this effect might by exaggerated by the interaction potentials used. This 
potential model does not give a protective van der Waals sphere to the charge on the hydrogen atom. In principle. 
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this is a fundamental difficulty, but in computer simulations, the heights of energetic barriers usually exclude the 
singularity. The development of interaction potentials for anion-water interactions nevertheless has to account for 
these problems. The strong interactions with the hydrogens "pull" the water closer and the first maxima of the 
ion-oxygen pair correlation function is already in the strongly repulsive region, reducing the effective ion radius. 

We have also studied the charging of sodium, potassium, calcium, fluoride, chloride, and bromide ions. The agree- 
ment with the available experimental data for solvation free energies is only qualitative, reproducing the trends with 
ionic size. The quantitative data are not in satisfactory agreement with the experimental results, even conceding quite 
substantial discrepancies between different compilations of the experimental data for certain ions. We observe typical 
errors of about 10 to 15 percent for the free energies of ionic solvation compared to the experimental data of Marcus.tS 
This clearly indicates the further need to develop quantitatively reliable descriptions of ion-water interactions. 

However, to allow for valid comparisons of data obtained from computer simulations with experimental results, it 
is crucial to eliminate systematic errors in the simulation methods. An important part of this study was devoted 
to analyzing the effect of finite system sizes on the free energy of charging. We could clearly establish that Ewald 
summation (and, similarly, the generalized reaction-field method) accounts for finite-size effects by adding a term that 
corrects for self-interactions. We showed that even for systems with only N — 16 water molecules it is possible to 
obtain accurate estimates of the solvation free energy of the sodium ion. For typical system sizes of a few hundred 
water molecules, these finite-size corrections are substantial in magnitude. Neglecting them yields results of little 
quantitative validity. 
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FIG. 1. Probability distributions P{e(j)) of the electrostatic energy e(j) at the site of a methane- like ion Me with charge q = e 

from Ewald summation with A*' = 256 (o, ), N = 128 (□, - • -), and GRF with TV = 256 water molecules (+, ), 

respectively. The lines are Gaussian distributions. Also shown are Gaussian distributions corrected for finite-size effects, which 
are peaked near e(f> — —550 kj mol~^; they agree closely in position and variance. 



12 



FIG. 2. The average electrostatic potential at the site of the negatively charged ion Me {q = — e) calculated from the pair 
correlations of a Monte Carlo simulation using Ewald summation and — 256 water molecules. The results of the integration 
using the GRF interaction with cutoff = L/2 and the bare Coulomb interaction 1/r are shown with long- and short-dashed 
lines, respectively. Finite-size corrections are added as constants. The Ewald-summation result is shown as a reference with a 
solid line. 



FIG. 3. The average electrostatic potential nic at the position of the methane-like Lennard- Jones particle Me as a function 
of its charge q. rric contains corrections for the finite system size. Results are shown from Monte Carlo simulations using Ewald 
summation with A'^ = 256 and A'^ = 128 (x) as well as GRF calculations with A'' — 256 water molecules (□). Statistical 
errors are smaller than the size of the symbols. Also included are linear fits to the data with q < and q > (solid lines). The 
fit to the tanh-weighted model of two Gaussian distributions (eq B2) is shown with a dashed line. 



FIG. 4. The fluctuation of the electrostatic potential /c at the position of a methane-like Lennard- Jones particle as a function 
of its charge q. fc contains corrections for the finite system size. Error bars indicate one estimated standard deviation of the 
data. For further details see Figure H. 



FIG. 5. The pair correlation functions gio (top panel) and gm (bottom panel) of the Me ion with water oxygen and hydrogen. 
The g{r) curves are shifted vertically according to the ionic charge by q/e, i.e., by 1 for q = e, 0.5 for q — 0.5e etc. The g{r) 
curves of Ewald summation and GRF simulations with A' = 256 water molecules are shown with solid and dashed lines, 
respectively. 



FIG. 6. The energetics of clusters of a fluoride ion and SPC water. Results are shown for the interaction energy Us of 
the fluoride ion with the water (o), as well as the mean (</>) (-I-) and variance (A(^^) (□) of the electrostatic potential at the 
ion position. The figure shows differences of these quantities with respect to the bulk values calculated from Monte Carlo 
simulation of an A'^ = 128 water-molecule system using Ewald summation: Aus — Ua — Us,ew , A{<f)) = {(j>) — ((^)ew, and 
A{A0^) = {A(/>^) — (A<^^)ew. The lines are fitted curves as explained in the text. Error bars indicate one standard deviation 
estimated from block averages. The standard deviations of the bulk and cluster data were added. 



TABLE I. Lennard- Jones parameters of the ion-water in- 
teractions. 



Ion e/(kJmol ^) cr/nm 

Na+ 0.200546 0.285000 

K+ 0.006070 0.452000 

Ca2+ 0.637972 0.317000 

F" 0.553830 0.305000 

or 0.537866 0.375000 

Br" 0.494464 0.383000 

Me 0.893228 0.344778 
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TABLE II. Results for the mean and the fluctuation of the 
potential 4> (with and without finite-size corrections) at the 
position of a methane-like Lennard- Jones particle Me carry- 
ing a charge q. "Coulomb" refers to the treatment of the 
electrostatic interactions (Ewald or GRF). A'^ is the number 
of water molecules. The mean and the fluctuation are listed as 
m = e{cj>) and / — f3e^{{(f>— {4>))^), both in units of kj mol"""^. 
The corrected values are nic = m + qe^ and fc ~ f — e^^. 
The statistical errors of m and / are estimated from block 
averages as approximately 4.0 and 30 kJ mol~^. 



N 


Coulomb 


g/e 


m 


/ 


rric 


fc 


256 


EW 


— 1.00 


670.3 


604 


869.9 


804 


128 


EW 


— 1.00 


618.1 


664 


869.2 


915 


256 


GRF 


-1.00 


520.5 


493 


869.1 


842 


128 


EW 


-0.75 


465.3 


730 


653.6 


981 


256 


EW 


-0.50 


324.0 


713 


423.7 


913 


128 


EW 


-0.50 


292.2 


698 


417.8 


950 


256 


GRF 


-0.50 


242.2 


568 


416.5 


917 


128 


EW 


-0.25 


141.3 


529 


204.0 


780 


256 


EW 


0.00 


38.0 


367 


38.0 


567 


128 


EW 


0.00 


37.3 


341 


37.3 


592 


256 


GRF 


0.00 


34.9 


239 


34.9 


587 


128 


EW 


0.25 


-40.7 


313 


-103.5 


564 


256 


EW 


0.50 


-143.1 


374 


-242.9 


573 


128 


EW 


0.50 


-118.1 


332 


-243.6 


583 


256 


GRF 


0.50 


-74.4 


232 


-248.7 


581 


128 


EW 


0.75 


-205.9 


354 


-394.2 


605 


256 


EW 


1.00 


-348.8 


450 


-548.3 


650 


128 


EW 


1.00 


-298.9 


389 


-550.0 


640 


256 


GRF 


1.00 


-202.2 


296 


-550.9 


645 



TABLE III. Free energy (in kJ moP^) of charging the 
methane-like Lennard- Jones particle Me from to ±e. The 
free energy was calculated from fltting to polynomials p„ of 
degree n and a tanh-weighted model of two Gaussian regimes 



(eqjl). 

Function A^'=^(0 ^ e) A^'=^(0 -e) 

'^2 -246 -454 

P4 -253 -431 

P6 -250 -431 

Pa -250 -431 

pio -250 -431 

tanh -250 -430 
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TABLE IV. Results for the mean rric and fluctuation /c of 
the potential (with finite-size corrections included) at the po- 
sition of sodium, potassium, calcium, fluoride, chloride, and 
bromide ions at different charge-states q. The data were cal- 
culated from Monte Carlo simulations using A'' = 128 water 
molecules and Ewald summation over 100 000 passes. The 
mean and the fluctuation are listed as rric = e{{(p) + qO 
fc = I3e'^{{(f)-{(f))f)-e'^i,, both in units of kj rnoPV The sta- 
tistical errors of rtic and fc are estimated from block averages 
to be approximately 4.0 and 30 kJ mol~^. 



Ion 


q/e 


rUc 


fc 


Na 


0.00 


39.0 


891 


Na 


0.50 


-395.6 


956 


Na 


1.00 


-885.1 


970 


K 


0.00 


38.6 


682 


K 


0.50 


-282.1 


690 


K 


1.00 


-662.6 


772 


Ca 


0.00 


41.0 


662 


Ca 


1.00 


-653.6 


789 


Ca 


2.00 


-1367.6 


667 


F 


0.00 


35.7 


718 


F 


-0.50 


587.6 


1381 


F 


-1.00 


1167.3 


961 


CI 


0.00 


36.2 


550 


CI 


-0.50 


378.2 


819 


CI 


-1.00 


794.1 


773 


Br 


0.00 


37.3 


545 


Br 


-0.50 


369.2 


758 


Br 


-1.00 


772.7 


773 
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TABLE V. Results for the free energy of charging the 
sodium cation from g = to e in SPC water in units of 
kj mol~^. ^j,"^ includes the finite-size corrections which are 
listed as Mscif. The free energies were calculated from poly- 
nomial fits to the derivative data of Table ^ (polynomials 
Pn of degree n). Also included are results of thermodynamic 
integration (TI). Linear charging paths from to e and from e 
to are denoted by upward (j) and downward (j) arrows, re- 
spectively. Ewald (EW) and generalized reaction-field (GRF) 
interactions were used for the charges. 



Method 


Coulomb 


A'" 






P2 


EW 


128 


-126 


-415 


Pa 


EW 


128 


—126 


-407 


Pa 


EW 


128 


—126 


-407 


Tlj 


EW 


256 


-100 


-404 


Tli 


EW 


256 


-100 


-406 


Tit 


EW 


128 


— 126 


-402 


Tli 


EW 


128 


-126 


-407 


TIT 


EW 


64 


-158 


-407 


Tli 


EW 


64 


-158 


-406 


TIT 


EW 


32 


-198 


-403 


Tli 


EW 


32 


-198 


-407 


TIT 


EW 


16 


-247 


-409 


Tli 


EW 


16 


-247 


-411 


TIT 


EW 


8 


-305 


-419 


Tli 


EW 


8 


-305 


-425 


TIT 


GRF 


128 


-219 


-401 


Tli 


GRF 


128 


-219 


-406 


TIT 


GRF 


64 


-276 


-408 


Tli 


GRF 


64 


-276 


-411 


TIT 


GRF 


32 


-346 


-419 


Tli 


GRF 


32 


-346 


-424 



TABLE VI. Results for the free energy ji^^ of charging the 
potassium, calcium, fiuoride, chloride, and bromide ions from 
q = to ±e, 2e in SPC water in units of kJ mol~^. ^"^^ 
includes finite-size corrections. Details as in Table M. 



Ion 


P2 


P4 


P6 


TIT 


Tli 


K+ 


-297 


-293 


-295 


-291 


-294 




-1317 


-1315 


-1316 


-1311 


-1327 


F" 


-594 


-590 


-590 


-590 


-594 


cr 


-401 


-392 


-392 


-389 


-394 


Br- 


-393 


-382 


-382 


-379 


-382 
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TABLE VII. Results for the excess chemical potential /i*^^ 
(in kj mol~^) of transferring an uncharged LJ particle from 
ideal gas into SPC water. The LJ parameters are those of 
Table M Errors are estimated from block averages. 



LJ-Particle 



Na 9.2(1) 

K 23.7(5) 

Ca 10.2(3) 

F 9.7(2) 

CI 21(3) 

Br 24(3) 

Me 10.2(9) 



TABLE VIII. Results for the calculated free energy of 
ionic hydration (in kJ moP^) compared with experimen- 
tal data. The experimental data were adjusted to give 
AG* = -1050 kJ mol"^ for H+, as used by Marcus.N Also 
included a*a computer simulatieii results by Straatsma and 



BerendseiH and Migliore et aZ.Ej 



Ion 


ex 


AG* " 


AG* 


AG* " 


ex d 

M 


ex C 

M 


ex f 


Na+ 


-398 


-365 


-371 


-372 


-508 


-431 


-459 


K+ 


-271 


-295 


-298 


-298 


-425 


-349 


-321 


Ce?+ 


-1306 


-1505 


-1553 




-1623 


-1394 




P- 


-580 


-465 


-394 


-441 


-497 


-421 


-418 


cr 


-371 


-340 


-277 


-324 


-315 


-239 


-237 


Br" 


-358 


-315 


-263 


-310 


-304 


-228 





^Experimental data of Marcus.S |— | 
'^Experimental data of Friedmapiand Krishnan.Eil 
'^Experimental data of Conway.E3 

''Computer simulation data of Straatsma and Berendsen 
calculated|— iising molecular dynamics of A'^ = 216 water 
molecules.LJ The results contain a Born-type correction ap- 
plied by the authors to their raw data. 

''Computer simulaticpidata of Straatsma and Berendsen with- 
out Born correction.^ 

'Computer simulation data of Migliore et al. calcji]|ated using 
molecular dynamics of A = 342 water molecules.Ej 
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Fig. 2 (Hummer, Pratt, and Garcia) 




Fig. 3 (Hummer, Pratt, and Garcia) 
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Fig. 4 (Hummer, Pratt, and Garcia) 



